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1.0  INTRODUCTION 

The  present  work  is  a result  of  efforts  directed  toward  providing  theoretical  support 
to  aerodynamic  testing  in  transonic  wind  tunnels.  The  corresponding  boundary-layer 
problems  usually  involve  compressible  turbulent  flows  with  both  favorable  and  adverse 
pressure  gradients  on  smooth  impermeable  adiabatic  walls.  The  following  requirements 
have  occurred  in  the  Propulsion  Wind  Tunnel  Facility  (PWT)  at  the  Arnold  Engineering 
Development  Center  (AEDC)  which  illustrate  the  need  to  solve  this  sort  of 
boundary-layer  problem  rapidly,  with  reasonable  accuracy,  and  with  minimal  preparation 
time:  (1)  support  of  tests  while  they  are  in  the  tunnel,  (2)  accounting  for  the 
boundary-layer  displacement  effect  in  solutions  requiring  viscous/inviscid  iterations,  and 
(3)  estimation  of  skin  friction  drag.  In  the  first  requirement,  quick  results  are  desired  to 
investigate  boundary-layer  problems  that  could  influence  the  remaining  test  plan.  In  the 
second  requirement,  transonic  inviscid  codes  can  be  complicated  and  require  large 
amounts  of  computer  time  and  storage;  hence  it  is  desirable  to  expend  minimum  time, 
effort,  and  computer  resources  on  the  boundary-layer  portion  of  the  problem.  For  the 
last  requirement,  it  is  usually  necessary  to  perform  computations  for  a large  number  of 
test  conditions  and  it  is  desirable  to  reduce  the  time  per  solution  to  conserve  computer 
time.  The  present  work  was  carried  out  to  yield  a turbulent  boundary-layer  calculation 
method  capable  of  providing  reasonably  accurate  solutions  in  short  periods  of  time  with 
minimal  demands  on  computer  resources. 

Turbulent  boundary-layer  computation  methods  can  be  divided  into  two  classes, 
integral  and  differential  techniques.  Integral  methods  require  the  solution  to  one,  or 
more,  ordinary  differential  equations,  but  differential  methods  require  the  solution  to  a 
system  of  partial  differential  equations.  As  a result,  development  of  differential  methods 
has  accelerated  in  recent  years  following  advances  in  computational  capabilities  (including 
both  computers  and  numerical  algorithms),  whereas  integral  methods  require  less  in  the 
way  of  computational  capability  and  have  been  under  investigation  for  many  years. 
Differential  methods,  in  conjunction  with  improved  turbulence  models  such  as  the  one- 
or  two-equation-of-turbulence  models,  are  usually  regarded  as  more  accurate  than  integral 
methods  for  general  applications  (Ref.  1).  With  regard  to  specific  applications,  however, 
the  results  of  the  Stanford  Conference  (Ref.  2)  on  incompressible  flow  indicate  that  some 
integral  techniques  provide  essentially  the  same  accuracy  as  differential  methods. 
Therefore,  an  integral  approach  is  considered  here  in  the  interest  of  simplicity  and  speed, 
with  the  additional  intent  of  developing  an  integral  method  to  the  point  of  being 
competitive  in  accuracy  with  differential  methods  for  transonic  turbulent  flow  over  an 
adiabatic  wall. 

The  approach  of  this  work  is  the  simultaneous  solution  of  the  momentum  integral 
equation  and  the  mean-flow  kinetic  energy  integral  equation  for  two-dimensional  planar 
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or  axisymmetric  flow.  The  momentum  integral  equation  is  usually  used  in  integral 
methods  (an  exception  to  this  is  White's  (Ref.  3)  method;  he  contends  that  his  technique 
is  entirely  different  from  von  Karman's  momentum  integral  equation).  Use  of  the 
mean-flow  kinetic  energy  integral  equation  is  not  uncommon  in  integral  techniques,  but  it 
is  not  universal.  A primary  advantage  of  using  this  latter  equation  was  pointed  out  by 
Nash  (Ref.  4)  to  be  that  it  is  extendible  to  compressible  flow  without  the  transformation 
limitations  associated  with  the  moment-of-momentum  integral  equation. 

In  addition  to  the  momentum  integral  equation  and  the  mean-flow  kinetic  energy 
integral  equation,  other  auxiliary  parameters  commonly  required  for  solution  of  the 
integral  equations  are  addressed: 

1 . The  wall  shear  stress 

2.  Velocity  profiles 

3.  The  shear-work  integral 

4.  The  velocity-temperature  relationship 

5.  Shape  factor  correlations 

The  methods  by  which  items  2 through  4 are  handled  in  the  present  work  are  different 
from  those  used  in  other  integral  techniques.  The  velocity  profile  is  a newly  developed 
profile  which  promises  increased  accuracy,  the  shear-work  integral  is  evaluated  directly 
rather  than  by  empirical  correlation,  ,and  the  temperature  (density)  distribution  is  based 
on  a new  relation  which  predicts  the  total-temperature  overshoot  near  the  outer  edge  of 
turbulent  adiabatic  wall  boundary  layers  with  nonunity  Prandtl  number  and  recovery 
factor. 

The  integral  method  is  developed  in  Section  2.0.  Numerical  solutions  are  compared 
with  experiments  and  finite  difference  boundary-layer  solutions  in  Section  3.0  to 
establish  the  capabilities  and  limitations  of  the  method.  These  solutions  are  all  for  an 
adiabatic  wall  as  this  is  the  situation  usually  encountered  in  transonic  flow  applications. 
The  resulting  computer  program  used  to  obtain  the  results  is  short  and  simple  to  use.  A 
listing  of  the  Fortran  code  is  presented  in  Appendix  A. 

2.0  ANALYSIS 

The  general  form  of  the  integral  boundary-layer  equations  for  x-momentum  and 
mean-flow  kinetic  energy  for  turbulent  compressible  flow  on  two-dimensional  planar  or 
axisymmetric  bodies  is  derived  in  this  section.  The  five  items,  which  amount  to  auxiliary 
relations,  mentioned  in  the  preceding  section  will  be  addressed  as  they  arise  in  the 
development  of  the  integral  method.  The  numerical  method  used  to  solve  the  equations  is 
a simple  one  used  by  Nash  (Ref.  4). 
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2.1  DERIVATION  OF  THE  INTEGRAL  EQUATIONS 


The  compressible  turbulent  boundary-layer  equations  of  continuity,  x-momentum, 
and  total  enthalpy  (energy),  are,  respectively, 


h (M  * £ (M  ' 0 


p u + pvis  = -in  + I d 
dx  dy  dx  rk  dy 


O) 

(2) 


pu 


dx 


u - pc 

dy 


(3) 


where  r = rw  + y cos  4/  for  external  flow,  r = rw  - y cos  0 for  internal  flow,  k = 0 for 
two-dimensional  planar  flow,  and  k = 1 for  axisymmetric  flow.  The  dependent  variables 
in  Eqs.  (1)  through  3 are  time-averaged  quantities  with  the  exception  of  the  term  pv, 
which  is  actually  the  product  of  the  time-averaged  p and  the  mass-weighted  averaged  v 
(see  Cebeci  and  Smith  (Ref.  5)  for  a discussion  of  the  approximations  leading  to  the 
boundary-layer  equations).  The  y-momentum  equation  is  dp/3y  = 0.  The  momentum  and 
mean-flow  kinetic  energy  integral  equations  will  be  derived  from  Eqs.  (1)  and  (2). 
Equation  (3)  is  the  basis  for  a relation  between  velocity  and  temperature  which  is 
required  for  solution  of  the  compressible  momentum  and  mean-flow  kinetic  energy 
integral  equations. 

The  momentum  and  mean-flow  kinetic  energy  integral  equations  are  derived  by 
multiplying  Eq.  (1)  by  um  + 1/(m  + 1),  Eq.  (2)  by  rkum  summing,  and  simplying  the 
results  to  obtain 


' _J_  l (rkpum  2)  + S (rkpvum~1)  = _rk  um  dp  + u">  d (r*  r) 

m — 1 dx  dy  dx  dy 


(4) 


where  r is  the  total  shear  stress 


r = p^-pu'v'  (5) 

Substitution  of  m = 0 will  eventually  result  in  the  momentum  integral  equation,  and  m = 
1 in  the  mean-flow  kinetic  energy  integral  equation.  Using  the  momentum  equation  at 
the  edge  of  the  boundary  layer 
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dp 

dx 


-P« 


du 

e 

dx 


(6) 


and  subtracting  the  term  l/(m  + 1)  9/9x  i*  pu  uem+1)  from  each  side  of  Eq.  (4),  one  can 
write  the  expression 


l 

m 1 


/um+1 

\u!r-+l 
v c 


( 


rkpvumTl 


)l 


- -L  / ^kP»«r‘V  u"  / 

dx  m + 1 dx  y f dy 


(rkr) 


(7) 


Expanding  the  middle  term  on  the  right-hand  side  of  Eq.  (7),  combining  part  of  this 
expansion  with  the  term  containing  due/dx,  and  using  Eq.  (1)  in  the  remaining  portion  of 
the  expansion  results  in  the  expression 


i 

m+  1 


Peu 


m+2 

e 


+ A.  /rk  p vum+1\  - um+1  A (rkpv) 
^ \ J e dy 


,m+l 
‘eJe 


(8) 


Integrating  Eq.  (8)  with  respect  to  y over  the  interval  (0  < y<  “)  yields  the  result 


The  no-slip  wall  boundary  condition  of  u = 0 at  y = 0 is  used  to  obtain  Eq.  (9).  In 
addition,  two  terms  resulting  from  the  integration  of  the  second  and  third  terms  on  the 
left-hand  side  of  Eq.  (8)  do  not  appear  in  Eq.  (9)  because  they  sum  to  zero  as  a 
consequence  of  pv  -*■  (pv)e  and  u-me  in  the  limit  of  y-*“. 


For  the  special  case  of  m = 0,  Eq.  (9)  can  be  written 


- rw  uePw  vw  = rw  rw  U0) 
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In  writing  Eq.  (10),  it  is  necessary  to  use  (i*t)-K)  as  y-*~.  To  show  that  (rkr)-K)  as  y-*-”, 
consider  the  following:  The  turbulence  model  used  subsequently  in  this  analysis  will 
postulate  the  shear  stress  in  the  outer  part  of  the  boundary  layer  to  be  proportional  to 
the  local  velocity  gradient.  Therefore,  because  r is  linear  in  y,  it  will  suffice  to  prove  that 
(y  3u/3y)  -*■  0 as  y-*~.  As  a consequence  of  the  velocity  .profile  developed  in  Ref.  6,  also 
to  be  used  subsequently  in  this  analysis,  the  term  3u/3y  behaves  as  exp  (-cyt> ) where  c > 
0 and  b > 0.  The  product  y exp(-cyb  )-*•  0 as  y-»~,  and  hence  (i*r)-K)  as  y-*“. 

For  the  special  case  of  m = 1,  Eq.  (9)  can  be  written 


DO 


0 


- dy 
P . ' 


DO 


(11) 


In  writing  Eq.  (11),  it  is  necessary  to  integrate  the  last  term  in  Eq.  (9)  by  parts  and  again 
use  (rk  r)-*0  as 


Interchanging  the  order  of  differentiation , and  integration  in  Eqs.  (10)  and  (11), 
these  equations  can  be  written  as 


where 


r p u 
w r e e 


2r*  p 
w r f 


cr  = r* 


2 p u2 
' e e 


s*  = 


e* 


e.  “I  0 * h.  £ 

e 

p V c . 

— r W W — 1 

p u 2 

r e e 

(12) 

u3  Q A + k ! - 

6 ' -e  d)t 

. ] Pvr  Vw  _ Cf  fl 

2 p ' u 2 

r e e 

(13) 

(14) 

o 

cos  dy 

(15) 

Pa  \(l  ± y cos 
Pt  Ue/  \ 

\ k 

dy 

(16) 

\ ^ 

cos  \!t  1 dy 

(17) 
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8‘*  = J JL  ^1  - ^1  - JL  cos  ^'k  dy  (18) 

D = f JL  g(u/ue)  (l  ± _y  COS  Ilf  )k  dy  (19) 

J T v <5y  \ rw  / 

o 

These  five  integral  quantities  are  denoted  by  the  tilde  (y)  in  order  to  distinguish  them 
from  true  defect  lengths.  For  example,  5*  defined  by  Eq.  (16)  is  the  true  mass-flow 
defect  (displacement)  length  for  two-dimensional  planar  flow  (k  = 0)  but  it  is  not  the 

true  mass-flow  defect  length  for  axisymmetric  flow  (k  = 1)  (see  Ref.  7).  However,  8* 

approaches  the  true  mass-flow  defect  length  for  axisymmetric  flow  when  the  boundary 
layer  is  thin  compared  to  the  local  body  radius.  This  is  true  not  only  because  y/rw  in  Eq. 
(16)  is  small  when  (l-pu/pe  ue)  is  largest,  but  also  because  6*/rw(a  term  not  appearing 
in  Eq.  (16)  but  that  does  appear  in  the  true  expression  for  the  mass-flow  defect  length, 
see,  e.g.,  Section  1.9  of  Ref.  5)  is  small  compared  to  unity.  It  is  important  to  make  this 
point  about  the  defect  lengths  because  if  Eq.  (16)  was  used  to  compute  the  displacement 
effect  on  an  axisymmetric  body  with  thick  boundary  layers  for  use  in  a viscous/inviscid 
iteration  calculation,  for  example,  the  computed  displacement  effect  would  be  incorrect. 
Although  the  lengths  defined  by  Eqs.  (15)  through  (18)  are  not  true  defect  lengths  for 
axisymmetric  flow,  Eqs.  (12)  and  (13)  are  exact  within  the  boundary-layer  approximation 
and  are  the  general  form  of  the  momentum  and  mean-flow  kinetic  energy  integral 
equations  for  compressible  flow  on  two-dimensional  planar  or  axisymmetric  bodies. 


The  analysis  is  hereafter  restricted  to  boundary  layers  that  are  thin  relative  to  the 
local  body  radius.  The  motivation  for  this  is  that  by  excluding  the  term  (y/rw)  cos  i//  in 
Eqs.  (15)  through  (18),  integrations  can  be  performed  once  and  for  all  in  order  to 
establish  correlations  for  certain  shape  factors  which  will  appear  in  the  analysis. 
Otherwise,  integration  of  Eqs.  (15)  through  (18)  has  to  be  performed  at  each  axial 
station  thereby  increasing  the  computational  time.  The  analysis  is  also  restricted  at  this 
point  to  an  impermeable  wall,  vw  = 0 It  would  be  a trivial  program  modification  to 
include  the  permeable  wall  capability,  but  the  associated  auxiliary  relations  were  not  in- 
vestigated as  to  their  suitability  for  use  in  permeable  wall  calculations.  For  thin  boundary 
layers  relative  to  the  local  body  radius,  and  for  impermeable  walls,  the  momentum  and 
mean-flow  kinetic  energy  integral  equations,  [Eqs.  (12)  and  (13)],  reduce  to 


2rk  p 

W r 


I JL  (rk  p u2  d)  - &L  du*  = Cf 

rfcp  a2  dx  V w Pe  e ) “ dT  T 

w r e e e 

_ A.  /Vk  p u3  + All  = c 

' dx  \ w Pe  e / n dx  : 


(20) 

(21) 
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where 


DO 

d = f Pa 

( l - a-  ^ dy 

(22) 

J Pe\ 

o 

V “e/ 

"7c- 

) * 

(23) 

oo 

6*  = f Pa 

J Pc  ue 

O 

(24) 

o 

(25) 

d.  r ^ 

<u/ue) 

dy 

(26) 

Equations  (20)  through  (26)  are  used  in  this  work,  where  k = 0 for  planar  flow  and  k = 
1 for  axisymmetric  flow.  (Actually,  in  the  computer  program,  planar  flow  is  handled  by 
inputting  rw  = constant.) 

2.2  WALL  SHEAR  STRESS 

An  expression  is  needed  for  the  local  skin  friction  coefficient  appearing  in  Eqs.  (20) 
and  (21)  in  terms  of  variables  appearing  in  these  equations.  Such  expressions  are 
frequently  derived  by  using  the  zero  pressure  gradient  (dp/dx  = 0)  momentum  equation 
to  develop  a transformation  relating  the  compressible  skin  friction  coefficient,  cf,  to  an 
incompressible  skin  friction  coefficient,  cf.  One  of  the  available  expressions  for  cf  can 
then  be  used  to  determine  cf.  In  the  present  work,  variable  pressure  gradient  (dp/dx  =£  0) 
as  well  as  zero  pressure  gradient  flows  are  of  interest  and  the  procedure  for  obtaining  an 
expression  for  cf  is  not  well  established.  Therefore,  the  rationale  behind  the  one 
developed  here  will  be  given;  the  success  of  the  resulting  expression  will  be  judged  by  the 
agreement  obtained  with  experimental  data. 

The  basis  for  the  present  expression  for  Cf  is  Coles'  law  of  corresponding  stations 
(Ref.  8)  which  is 


cf 


R*e 


Pe  Pe 
PWPW 


cf 


(27) 


where  the  superscript  bar  (”)  denotes  low-speed  flow.  This  relation  was  used  by  Lewis, 
Kubota,  and  Webb  (Ref.  9),  along  with  the  approximation,  n ~ T,  to  successfully 
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correlate  high-speed  (up  to  Me  - 8.18)  turbulent  boundary-layer  velocity  profiles.  Using 
this  approximation,  i.e.,  n ~ T,  Eq.  (27)  reduces  to 

Cf  Re0  = (28) 

It  was  pointed  out  in  Refs.  8 and  9 that,  whereas  Eq.  (27)  specifies  the  relation  for 
corresponding  stations,  it  does  not  specify  the  relation  between  cf  and  cf,  which  must 
come  from  hypothesis.  In  this  work,  recourse  is  made  to  the  work  of  Winter  and  Gaudet 
(Ref.  10),  who  established  the  relation 


for  adiabatic  walls,  dp/dx  - 0,  and  air.  It  is  hypothesized  that  Eq.  (29)  is  acceptable  for 
dp/dx  9^0  if  an  appropriate  expression  is  used  for  cf  for  dp/dx  =£  0.  The  expression 
used  for  cf  is 


/,  — .1.74+0.3 1H 

Uog10  Rep) 

which  is  Eq.  (6-179)  in  White's  book  (Ref.  3).  The  incompressible  momentum  thickness 
Reynolds  number,  Ree,  used  in  Eq.  (30)  is  obtained  using  Eqs.  (28),  (29),  and  Re„ 
which  is  calculated  as  part  of  the  solution.  The  incompressible  shape  factor,  H,  used  in 
Eq.  (30)  is  obtained  by  a subsequently  developed  correlation  of  H as  a function  of  H6  * 
(which  is  also  calculated  as  part  of  the  solution)  and  Me. 

2.3  VELOCITY  PROFILE 

The  integrals  of  Eqs.  (22)  to  (25)  are  usually  evaluated  by  use  of  power-law  or 
law-of-the-wake  (Ref.  11)  velocity  profiles.  The  main  feature  distinguishing  the  current 
method  from  previous  integral  methods  and  the  main  contribution  to  the  quality  of  the 
results  is  the  use  of  a newly  developed  velocity  profile.  This  new  profile  is  given  by  a 
closed-form  analytical  expression  valid  over  the  entire  thickness  of  the  turbulent 
boundary  layer.  Ref.  6.  Shown  in  step  8 of  Table  1,  the  expression  is  the  sum  of  two 
transcendental  functions,  one  expressed  in  terms  of  the  inner  variable,  y+,  and  one 
expressed  in  terms  of  the  outer  variable,  y/0.  Table  1 was  taken  from  Ref.  6 and  is  a 
summary  of  the  procedure  for  calculating  turbulent  boundary-layer  velocity  profiles  by 
this  method.  This  expression  was  shown  in  Ref.  6 to  provide  a reasonably  accurate  means 
of  correlating  various  velocity  profiles  that  included  boundary  layers  near  separation, 
reattached  boundary  layers,  nonequilibrium  boundary  layers  (boundary  layers  with 
nonzero  streamwise  variation  of  0),  and  compressible  boundary  layers.  In  the  latter  case, 
the  incompressible  profile  given  in  Table  1 is  related  to  the  corresponding  compressible 
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profile  by  u/ue  = u/ue  and  y/0  = (y/0)(0/6>)  where  0/0  is  a correlation  developed  in 
Subsection  2.6.  It  is  convenient  that  the  expression  for  the  profile  is  a function  of  the 
local  incompressible  parameters:  skin  friction  coefficient,  shape  factor,  and  Reynolds 
number  based  on  momentum  thickness,  because  these  parameters  are  a natural 
consequence  of  the  equations  being  solved.  This  velocity  profile  expression  was  developed 

Table  1.  Summary  of  Procedure  for  Computation  of  Turbulent 
Boundary- Layer  Velocity  Distributions 


Requirement 


H,  u^,  and  Re^  are  Inputs  (use  Incompressible  values) 


Compute  — (2)  = 


1.723  e-0-6  H^1  + 


Compute  — (5)  - 0.87  + 0 ,08e-2 ' 6(B_1 ,95^ 


Compute  g{2) 


(2) i-—  tan-1 

u 0.09  u+ 

e e 


0*18  \T 


— (5) 1 tan-1 


Compute  g(5) 


0.09  ul 


(0.18  Re A 
> ^ / 

ts 


LfflatZ.I 


H_  n\  at  y - 5 


ue  \A 


gA  at  y - 2 


g W at  L - 5 


In  ftanh-1  [g2(2)' 

Compute  b - ^tanh .g 

Compute  a = tanh  1 [g2(2>] 

2b 


tanh  1z  - i In 


tanh  z 


H H_ 

u_  - 


e2z  - 1 
e2z  + 1 


8-  3759  tan-^o.oe  y+)  - (K  - rig)  tanh1/2  ag 


ue  K 


for  two-dimensional  planar  flows,  but  it  is  used  herein  for  axisymmetric  flows  without 
modification,  because  the  assumption  has  been  made  that  the  boundary  layer  is  thin 
relative  to  the  local  body  radius. 
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Because  the  new  profile  reduces  to  an  identity  at  the  wall,  it  is  not  possible  to 
calculate  the  wall  shear  stress,  rw  = (ju9u/9y)w , by  taking  the  slope  of  the  velocity  profile. 
Instead,  the  wall  shear  is  obtained  from  the  skin  friction  coefficient  developed  in 
Subsection  2.2.  The  shear  work  integral,  Eq.  (26),  is  usually  obtained  by  correlation  of  , 
experimental  data,  but  the  new  velocity  profile  is  considered  sufficiently  accurate  to 
permit  evaluation  of  Eq.  (26)  by  use  of  9u/9y  obtained  from  the  profile  expression  and 
an  appropriate  turbulence  model,  as  described  in  Subsection  2.4. 

2.4  SHEAR-WORK  INTEGRAL 


The  term  referred  to  here  as  the  shear-work  integral  (also  frequently  referred  to  in 
the  literature  as  the  dissipation  integral)  is  denoted  by  D and  given  by  Eq.  (26).  The 
velocity  profile  given  in  Table  1 is  used  to  compute  the  velocity  derivative,  and  in  this 
regard  it  is  convenient  to  express  Eq.  (26)  as 


D = 


<HU  ■ ue) 

<5(y  / 9 ) 


cl  (y/0  ) 


(31) 


inasmuch  as  y/0  is  one  of  the  independent  variables  in  the  velocity  expression.  The 
problem,  therefore,  is  to  specify  the  shear-stress  distribution,  r/rw,  across  the  boundary 
layer  in  order  to  carry  out  the  integration  identified  by  Eq.  (31).  In  keeping  with  the 
idea  of  maintaining  simplicity,  the  shear-stress  distribution  is  specified  in  three  regions 
across  the  boundary  layer,  and  common  turbulence  models  for  each  region  are  used. 


The  first  region  considered  is  that  near  the  wall,  defined  by  0 < y+  < 100.  In 
deriving  the  inner  solution  part  of  the  velocity  profile  in  Table  1 (the  inner  solution  part 
is  tan-l  (0.09  y+)/0.09),  it  was  assumed  that  the  total  shear  stress  r was  constant  and 
equal  to  rw . Based  on  the  agreement  between  experimental  data  and  the  velocity  profile 
expression  resulting  from  this  assumption  (see  Ref.  6),  it  appears  that  t = tw  for  0 < y+ 

< 100  is  a plausible  approximation  for  use  in  Eq.  (31).  In  fact,  in  Fig.  7 of  Ref.  6 an 
extreme  case  of  an  experimental  velocity  profile  having  a severe  adverse  pressure  gradient 
and  no  logarithmic  region  was  well  represented  by  the  velocity  expression  in  Table  1 , and 
for  this  case  the  inner  solution,  tan-l  (0.09  y+)/0.09,  was  the  dominant  term  for  0 < y+ 

< 100.  The  assumption  t = rw  allows  Eq.  (31)  to  be  integrated  analytically  over  0 < y+ 

< 100  (i.e.,  out  to  the  value  of  y/0  where  y4  = 100  as  given  by  y/0  = 100  ue+/RJa,  see 
Table  1)  and  numerical  integration  is  avoided. 


The  next  two  regions  axe  handled  by  Boussinesq's  (Ref.  12)  eddy-viscosity  concept, 
where  the  eddy  viscosity  e is  defined  by 


/ * 

— n v = 


(32) 
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where  e will  be  either  em  or  e0,  depending  on  the  region.  For  the  first  region  outside  y+ 
= 100  (middle  region),  Prandtl's  mixing-length  theory  (Ref.  12)  is  used  where 


and  the  mixing  length  is 


<?u 

dy 


£=  0.41  y 


(33) 


(34) 


For  the  outer  region  a constant  eddy  viscosity  of  the  form  suggested  by  Clauser  (Ref.  13) 
is  used.  The  specific  expression  for  e0  in  the  outer  region  is 


f0  = 0.0168  ue  8* 


where  5*  is  the  incompressible  displacement  thickness  defined  as 

M"  Hb 


(35) 


(36) 


The  point  in  the  boundary  layer  separating  the  middle  region  from  the  outer  region  at  a 
particular  x location  is  that  point  where  em  = e0.  The  shear  stress  for  y+  > 100  primarily 
consists  of  the  turbulent  shear  stress  as  defined  by 


(37) 


and  the  shear  stress  distribution  for  use  in  Eq.  (31)  is  approximated  as 


p (0.41  y) 2 

du 

du 

dy 

dy 

-2_  _P  / 0.41  X\2 

ci  \ 6 ) 


Xf) 


Xi) 


(38) 


for  the  middle  region,  and 


5ll 

0.0i68p u 5* -5—  \ 

JL  = 3y  = 0.0168  J_  £-  H £ (u/ue ) 


cf  Pe 


d(y/6  ) 


(39) 


for  the  outer  region. 


Equation  (31)  is  analytically  integrated  from  y+  = 0 to  y+  = 100  and  then 
numerically  integrated  in  the  middle  and  outer  regions,  using  Eqs.  (38)  and  (39)  with 
p/pe  being  evaluated  as  discussed  in  Subsection  2.5. 

2.5  VELOCITY-TEMPERATURE  RELATIONS 

To  establish  the  necessary  shape  factor  correlations  for  use  in  the  solution  to  Eqs. 
(20)  and  (21)  it  is  necessary  to  use  a relation  between  velocity  and  temperature  in  order 
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to  integrate  Eqs.  (22)  to  (25)  (recall  that  plpt  = Te/T  as  a consequence  of  dp/dy  = 0). 
The  relationship  used  is  the  analytical  solution  developed  in  Ref.  14  for  nonunity  Prandtl 
number.  Because  the  interest  here  is  in  turbulent  boundary  layers  on  adiabatic  walls,  only 
the  adiabatic  wall  solution  in  Ref.  14  is  used.  The  advantage  of  the  method  developed  in 
Ref.  14  over  the  more  commonly  used  method  of  Crocco  (Ref.  15)  is  that  the  results 
obtained  in  Ref.  14  predict  the  overshoot  in  total  temperature  near  the  outer  edge -of  a 
turbulent  boundary  layer  on  an  adiabatic  wall  [see,  e.g.,  van  Driest  (Ref.  16)],  which  the 
Crocco  results  cannot  do.  Moreover,  the  Crocco  solution  is  actually  contained  as  a special 
case  of  the  solution  given  in  Ref.  14  for  unity  Prandtl  numbers.  The  quality  of  agreement 
between  the  results  of  Ref.  14  and  experimental  data  is  considered  good.  The 
comparisons  from  which  this  conclusion  was  drawn  are  given  in  Refs.  14  and  17. 

The  density  ratio  in  the  expressions  for  the  shear-stress  distributions,  Eqs.  (38)  and 

(39),  is,  however,  specified  by  the  modified  Crocco  relation  for  an  adiabatic  wall  (see, 
e.g.,  Ref.  3) 

pe  t p = T/'Te  = 1 + r J{y  - 1)  (l  - u2  / Ug  )/2 

The  reason  for  this  choice  rather  than  using  Ref.  14  is  that  numerical  integration  of  Eq. 
(31)  is  carried  out  at  each  axial  location  (whereas  integrations  for  the  shape  factor 
correlations  are  carried  out  only  one  time)  and  the  modified  Crocco  relation  is  simpler 
than  the  solution  given  in  Ref.  14,  requiring  less  computer  time.  Moreover,  the  modified 
Crocco  relation  is  in  good  agreement  with  the  solution  given  in  Ref.  14  in  the  region  near 
the  wall  for  adiabatic  flow  and  this  region  is  the  dominant  portion  of  the  integrand  of 
Eq.  (31). 

2.6  SHAPE  FACTOR  CORRELATIONS 

There  are  four  length  scales  in  Eqs.  (20)  and  (21),  0,  5*,  0*,  and  5**.  It  is 
convenient  and  computationally  expedient  to  correlate  these  lengths  by  analytical 
expressions  of  other  variables  to  avoid  evaluation  of  the  integrations  identified  by  Eqs. 
(22)  through  (25)  at  each  x location  (and  possibly  several  times  at  each  x location  to 
obtain  convergence).  The  approach  is  to  define  the  following  shape  factors 


(40) 

"e* ' f 

(41) 

H«"  ' T 

(42) 
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and  then  correlate  these  shape  factors  in  terms  of  Me  and  H where 

u = (43) 

0 

In  addition,  because  y is  scaled  by  6 as  shown  in  Table  1,  and  because  6 is  obtained  from 
the  solution  to  Eqs.  (20)  and  (21),  a correlation  for  did  is  also  established. 

The  velocity  profile  given  in  Table  1 was  used  to  carry  out  the  numerical 
integrations  (using  Simpson's  rule) ^necessary  to  establish  the  correlations.  This  profile 
depends  on  cf  and_Ree  as  well  as  H.  Using  Eq.  (30)  for_cf  = cf  (Re*,  H),  the  variables 
are  reduced  to  Re0  and  H only.  The  influence  of  Ree  on  the  numerical  results  is 
extremely  small  compared  to  that  of  H.  For  example,  the  difference  between  solutions 
with  Re0  - 5,000  and  Re9  = 100,000  is  less  than  one  percent,  which  is  within  the 
accuracy  of  the  analytical  expressions_derived  for  the  correlations.  All  numerical  results 
given  here  for  the  correlations  are  for  Ree  = 50,000. 

The  numerical  results  for  Hs*,  He*,  H6**,  and  dfd  are  represented  in  Figs.  1 
through  4 as  open  symbols.  The  results  of  the  analytical  expressions  derived  to  represent 
these  numerical  results  are  given  in  Figs.  1 through  4 as  solid  lines.  The  analytical 
expressions  derived  by  inspection  are 


Z , /r  %’ 

Id  ^ 
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and 


1I5*  = H (l  + 0.113  mJ  ) + 0.290  Me2 


0.92  M2 

e 

7.09  4-  M2 
e 


(44) 

(45) 

(46) 

(47) 


Although  the  constants  differ,  the  general  form  of  Eqs.  (44)  and  (46)  was  used  by  Nash 
(Ref.  4)  for  these  shape  factor  correlations.  The  form  of  Eq.  (45)  was  not  used  by  Nash, 
however.  The  advantage  of  Eq.  (45)  over  the  form  used  by  Nash  is  that  Eq.  (45)  can  be 
inverted  to  obtain  H as  a function  of  He*  and  Me,  which  is  useful  in  the  numerical 
solution  of  Eqs.  (20)- and  (21).  The  correlation  given  by  Eq.  (47)  is  in  good  agreement 
with  that  established  by  Winter  and  Gaudet  (Ref.  10),  although  their  correlation  was  for 

i* 


\ 
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Open  Symbols  Represent  the  Numerical  Results  for 
Etj.  (40)  Using  the  Velocity  Profile  Given  in  Table  1 
and  the  Velocity-Temperature  Relation  of  Ref.  14 

Solid  Curves  Represent  the  Analytical  Correlation  of 
H6*  = H (1  + a 113  m|)  + 0. 290  [Eq.  (44)1 


Figure  1.  Correlation  of  the  shape  factor  Hg.  for  adiabatic  flow. 


Open  Symbols  Represent  the  Numerical  Results  for 
Eq.  (41)  Using  the  Velocity  Profile  Given  in  Table  1 
and  the  Velocity-Temperature  Relation  of  Ref.  14 
Solid  Curves  Represent  the  Analytical  Correlation  of 
, 2 - 0l  54  tanh  ( 1. 1 (H  - 1)]  + 0.028m£ 

IEl“511 


Figure  2.  Correlation  of  the  shape  factor  Hfl*  for  adiabatic  flow. 
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Open  Symbols  Represent  the  Numerical  Results  for 
Eq.  (42)  Using  the  Velocity  Profile  Given  in  Table  1 
and  the  Velocity-Temperature  Relation  of  Ref.  14 

Solid  Curves  Represent  the  Analytical  Correlation  of 
• + 0. 25lj  A/l|  [Eq.  (46)1 


Figure  3. 


Correlation  of  the  shape  factor  Hg*  • 


for  adiabatic  flow. 


Open  Symbols  Represent  the  Numerical  Results  for 
0/9  Using  the  Velocity  Profile  Given  in  Table  1 and 
the  Velocity-Temperature  Relation  of  Ref.  14 
Solid  Curves  Represent  the  Analytical  Correlation  of 
a 0. 92  A/I? 

4-  - 1 V tanh  1 1.49  (H  - 0. 9)1  (Eq.  147)1 

0 7.09+1 Vl| 


Figure  4.  Correlation  of  the  compressible  to  incompressible 
momentum  thickness  ratio  for  adiabatic  flow. 

zero  pressure  gradient  flows  and  did  not  contain  the  parameter  H.  The  correlations  given 
by  Eqs.  (44)  through  (47)  and  compared  with  numerical  results  in  Figs.  1 through  4 are 
for  Me  < 3 and  a constant  recovery  factor  of  0.88. 

2.7  SUMMARY  OF  THE  EQUATIONS  AND  METHOD  OF  SOLUTION 

Using  the  definitions  of  the  shape  factors  defined  by  Eqs.  (40)  through  (42),  Eqs. 
(20)  and  (21)  can  be  rewritten  as 
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^3  * ♦■w  £-?» 


w r e e 


(48) 

(49) 


where  Hs  * , H„  * , and  H6  * • are  given  by  Eqs.  (44)  through  (46)  and  are  functions  of  H 
and  Me  as  indicated.  The  skin  friction  coefficient,  cf,  used  in  the  solution  of  Eqs.  (48) 
and  (49)  is  given  by  Eqs.  (28)  through  (30)  and  can  be  written  as 


cf  = J 
1 F 


0.3e 


-1.33H 


C 0*8  10  (»VFc>]  1-7 


74-0.3  1H 


(50) 


where  Fc  = (1  + M^/5)1 1,2 . The  skin  friction  coefficient  is,  therefore,  a function  of  H, 
Me,  and  Re0  where  Re0  = peue  0//ie  The  shear-work  integral,  D,  is  defined  by  Eq.  (31) 

where  t/t  w is  unity  for  0 < y+  < 100,  and  given  by  Eqs.  (38)  and  (39)  for  the  middle 
and  outer  regions  of  the  boundary  layer.  Therefore,  D depends  on  cf,  H,  6,  and  p/pe 
(where  pe/p  is  given  by  the  modified  Crocco  relation.  Subsection  2.5).  All  of  these 
quantities  can  be  related  to  H,  Me,  and  Ree  (or  6 in  place  of  Re0  because  peue/pe  is 
known  from  the  boundary-layer  edge  conditions).  Equations  (48)  and  (49)  are,  therefore, 
two  equations  for  the  two^  unknowns  H and  6 . 


The  method  used  to  solve  Eqs.  (48)  and  (49)  is  the  predictor-corrector  method  used 
by  Nash  (Ref.  4).  Various  schemes  were  tried  for  determining  step  size  and  for  most 
problems  the  step  size  was  not  particularly  important.  A constant  value  for  the  step  size 
was  used  for  the  solutions  presented  in  this  report;  however,  because  the  code  is  simple, 
it  is  a trivial  task  for  the  user  to  modify  this  portion  of  the  program  as  appropriate. 
Initial  values  of  6 and  H (from  which  H is  obtained  in  the  program  using  Eq.  (44))  must 
be  input  as  well  as  a reference  Mach  number,  Reynolds  number,  and  Me  and  rw  as 
functions  of  x.  A listing  of  the  boundary-layer  code  is  given  in  Appendix  A. 

3.0  RESULTS 


Because  a turbulent  boundary-layer  calculation  method  can  be  made  to  fit  a 
particular  set  of  experimental  data,  it  is  of  interest  to  consider  several  sets  of  data  for 
various  conditions.  The  present  method  was  applied  to  a large  number  of  experiments 
and  examples  are  presented  in  this  section.  The  examples  presented  are  those  of 
Schubauer  and  Klebanoff  (Ref.  18),  Stratford  (Ref.  18),  Cook,  McDonald,  and  Firmin 


s 
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(Ref.  19),  Pasiuk,  Hastings,  and  Chatham  (Ref.  20),  and  Winter,  Rotta,  and  Smith  (Ref. 
21).  These  examples  represent  a fairly  broad  range  of  flow  conditions  that  include 
favorable,  mild,  and  severe  adverse  pressure  gradients. 

The  experiment  of  Schubauer  and  Klebanoff  (Ref.  18)  was  selected  because  it  is  a 
realistic  airfoil-type  flow  commonly  encountered  and  also  numerous  other  calculation 
techniques  were  compared  to  this  experiment  in  Ref.  2,  from  which  indirect  comparisons 
with  the  present  method  can  be  inferred.  Comparisons  between  calculated  and  measured 
values  of  cf,  Hg*.  and  Re0  are  given  in  Fig.  5.  The  agreement  is  reasonably  good 
upstream  of  about  x = 22.5  ft.  This  is  the  same  sort  of  agreement  obtained  by  nearly  all 
the  27  calculation  techniques  (both  integral  and  finite  difference)  shown  in  Ref.  2. 


Figure  5.  Present  calculations  and  incompressible  experimental 
measurements  of  Schubauer  and  Klebanoff  (IDENT  = 
2100,  Ref.  18). 
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However,  those  calculations  shown  in  Ref.  2 which  took  three-dimensional  flow  effects 
into  account  were  in  better  agreement  with  these  experimental  data  than  the  present 
calculations.  Three-dimensional  effects  were  not  considered  in  the  present  calculations. 
Comparisons  of  calculated  and  measured  velocity  profiles  are  shown  in  Fig.  6 for 
favorable,  nearly  zero,  and  severe  adverse  pressure  gradient  portions  of  the  flow  as 
indicated  by  the  experimental  values  of  /3exp.  This  agreement  is  considered  good,  even  for 
the  severe  adverse  pressure  gradient  at  x = 22.0  ft. 


Figure  6.  Favorable  and  adverse  pressure  gradient  incompressible  velocity 
profiles  of  Schubauer  and  Klebanoff  (IDENT  = 2100,  Ref.  18). 

In  contrast  to  the  experiment  of  Schubauer  and  Klebanoff  (Ref.  18)  where  nearly 
all  27  calculation  techniques  (Ref.  2)  were  able  to  compute  the  flow  with  reasonable 
success,  the  experiment  of  Stratford  (Ref.  18)  had  comparisons  presented  from  only 
about  half  of  these  calculation  techniques  and  of  these  only  two  succeeded  in  calculating 
reasonably  good  results  for  cf,  Hfi  * , and  Re0  (Ref.  2).  Comparisons  between  the  present 
calculations  and  the  experiment  of  Stratford  (Ref.  18)  are  given  on  Figs.  7 and  8.  Of 
particular  interest  in  these  comparisons  is  the  velocity  profile  in  Fig.  8 for  /3cxp  = 78.988 
which  appears  to  be  extremely  close  to  separation  and  yet  the  experimental  velocity 
distribution  is  reasonably  well  calculated. 

Boundary-layer  measurements  on  an  RAE  2822  airfoil  in  transonic  flow  have 
recently  been  made  by  Cook,  McDonald,  and  Firmin  (Ref.  19).  Three  sets  of  data  from 
Ref.  19  are  considered,  and  the  pressure  distributions  are  given  in  Fig.  9.  Case  2 
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X,  ft 

Figure  7.  Present  calculations  and  incompressible  experimental 
measurements  of  Stratford  (I DENT  = 5300,  Ref.  18). 

represents  a reasonably  mild  favorable  and  adverse  pressure  gradient  flow  where  the  flow 
is  subcritical  and  the  airfoil  is  at  a = -2.18  deg.  Cases  9 and  12  represent  commonly 
encountered  supercritical  flows  with  severe  adverse  pressure  gradients  caused  by  shocks 
on  the  airfoil  (a  = 3.19  deg).  Comparisons  between  calculated  and  measured  values  of  cf, 
S*/c,  9/c,  and  H6*  are  given  in  Fig.  10  for  Case  2.  These  quantities  vary  smoothly  along 
the  airfoil  and  the  agreement  between  calculations  and  measurements  is  good. 
Comparisons  of  the  same  quantities  are  made  for  Case  9 in  Fig.  11.  Good  agreement  is 
obtained  up  to  the  shock  and  reasonable  agreement  is  obtained  downstream  of  the  shock. 
Case  12  was  for  the  same  conditions  as  Case  9 except  the  chord  Reynolds  number  was 
reduced  from  6.5  x 106  to  2.7  x 106.  Comparisons  of  calculated  and  measured  data  are 
presented  in  Fig.  12  for  Case  12.  The  agreement  in  Fig.  12  is  good  except  for  the  integral 
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quantities  at  x/c  = 0.9.  This  experiment  (Ref.  19)  provides  the  opportunity  to  investigate 
the  shock/boundary-layer  interaction  region  in  some  detail.  The  three  velocity  profiles 
nearest  the  interaction  (which  are  located  just  upstream,  essentially  at,  and  just 
downstream  of  the  interaction)  for  Cases  9 and  12  are  compared  with  the  calculations  in 
Figs.  13  and  14.  The  agreement  in  Fig.  13  for  the  higher  Re^g  flow  with  an  abrupt 
change  in  boundary-layer  properties  (see,  e.g.,  the  cf  distribution  in  Fig.  11  which 
illustrates  this  abrupt  change)  is  considered  reasonably  good.  The  agreement  in  Fig.  14  is 
considered  good  with  only  a slight  discrepancy  between  theory  and  experiment  for  the 
velocity  profile  in  the  interaction  region  at  x/c  = 0.498. 


o Experimental  Data,  IDENT  = 5300,  Ref.  18 
x- 4. 103  ft,  Pexp=31.320 


□ Experimental  Data,  IDENT -5300,  Ref.  18 
x - 6. 236,  Pexp  - 78. 988 


U/Ug 


Figure  8.  Severe  adverse  pressure 
gradient  incompressible 
velocity  profiles  of 
Stratford  {IDENT  = 
5300,  Ref.  18). 


Sym  ^co  rcem,  c*1q6  q,  deg  Case 

O 0.676  5.7  -2.18  2 

o 0.730  6.5  3.19  9 


0 02  0.4  0.6  0.8  1.0 

x/c 


Figure  9.  Upper  surface  pressure 
distributions  on  the 
RAE  2822  airfoil  for 
Cases  2,  9,  and  12  of 
Ref.  19. 


Measurements  by  Pasiuk,  Hastings,  and  Chatham  (Ref.  20)  were  made  for  a 
compressible  turbulent  adiabatic  wall  boundary  layer  with  a favorable  pressure  gradient. 
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Boundary-layer  displacement  and  momentum  thickness  measurements  are  compared  to 
present  calculations  in  Fig.  1 5 and  the  agreement  is  considered  good.  The  boundary-layer 
edge  Mach  number  varied  from  about  1.5  to  3.0  for  the  data  shown  in  Fig.  15. 
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Figure  10.  Upper  surface  boundary-layer 
properties  of  the  RAE  2822 
airfoil  for  = 0.676,  Re^c 
5.7  x 106,  and  a = -2.18  deg 
(Case  2). 


Figure  11.  Upper  surface  boundary-layer 
properties  of  the  RAE  2822 
airfoil  for  M,,,  = 0.730, 

Re^  c = 6.5  x 106,  and  a = 
3.19  deg  (Case  9). 


An  experimental  investigation  was  performed  by  Winter,  Rotta,  and  Smith  (Ref.  21) 
on  a waisted  body  of  revolution  in  both  subsonic  and  supersonic  flow.  Two  sets  of  data 
from  Ref.  21  are  considered,  M = 0.597,  Re  { = 9.98  x 106  and  M = 1.398,  Re  p = 
10.08  x 106.  The  pressure  distributions  for  these  two  cases  and  the  body  geometry  are 
shown  in  Fig.  16.  Boundary-layer  properties  of  cf,  6*,  0,  and  Hs  * are  compared  with 
present  calculations  in  Fig.  17  for  = 0.597  and  in  Fig.  18  for  = 1.398.  The 
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calculated  cf  is  a little  larger  than  the  experiment  for  the  = 0.597  data,  but  the 
overall  agreement  in  Fig.  17  is  considered  good.  The  agreement  with  the  = 1.398  data 
in  Fig.  18  is  good  up  to  x = 42  in.,  but  not  as  good  over  the  remaining  aft  end  of  the 
body.  In  this  connection,  a compatibility  check  of  the  experimental  data  was  performed 
by  Winter,  Rotta,  and  Smith  (Ref.  21)  by  using  measured  values  of  cf  and  H6*  in  the 
momentum  equation.  The  conclusion  drawn  in  Ref.  21  from  this  exercise  was  that  for  the 
Ho  = 1 -398  data,  the  skin  friction  measurements  appear  to  be  too  large  for  x > 42  in. 
This  was  attributed  in  Ref.  21  to  the  technique  used  to  measure  skin  friction  (the  razor 
blade  technique)  which  had  not  been  completely  established  for  these  types  of  flow. 


x/c 


Figure  12.  Upper  surface  boundary-layer 
properties  of  the  RAE  2822 
airfoil  for  = 0.730,  ReM  c 
= 2.7  x 106,  and  a = 3.19 
deg  (Case  12). 
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Figure  13.  Velocity  profiles  in  the  shock  region  of  the  RAE  2822 
airfoil  for  M,  = 0.730,  Re^,.  = 6.5  x 106,  and 
a = 3.19  deg  (Case  9). 


u/ue 

Figure  14.  Velocity  profiles  in  the  shock  region  of  the  RAE  2822 
airfoil  for  M„  = 0.730,  Re  . = 2.7  x 106,  and 
a = 3.19  deg  (Case  12). 
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X,  ft 

Figure  15.  Boundary-layer  displacement  and  momentum 
thickness  distributions  for  a compressible 
adiabatic  wall  boundary  layer  with  a 
favorable  pressure  gradient. 

Velocity  profiles  for  these  two  sets  of  data  from  Ref.  21  are  compared  with  present 
calculations  in  Figs.  19  and  20.  The  agreement  is  good  except  for  the  last  two 
downstream  profiles  in  Fig.  20  which  is  poor.  This  is  the  same  region  of  discrepancy  as  in 
Fig.  18.  Also  shown  in  Figs.  19  and  20  are  the  ratios  of  S*/rw  which  must  be  small  (i.e., 
6*/rw  «1)  for  transverse  curvature  to  be  unimportant  (recall  that  this  method  is 
axisymmetric  but  does  not  account  for  transverse  curvature).  Note  that  regions  of  any 
significant  discrepancies  in  Figs.  19  and  20  (excluding  the  last  two  downstream  profiles  in 
Fig.  20)  occur  near  the  outer  edge  of  the  boundary  layers.  It  was  pointed  out  in  Ref.  6, 
where  the  velocity  profile  used  here  was  developed,  that  the  correlation  for  u/ue  at  y/0  = 
5 should  be  modified  for  axisymmetric  flow.  Therefore,  the  strictly  two-dimensional 
planar  correlation  at  y/8  = 5 is  identified  as  the  source  of  this  discrepancy  in  Figs.  19 
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and  20.  It  should  be  pointed  out,  however,  that  the  portion  of  the  boundary  layer  for 
y Id  < 5 is  the  important  part  for  the  calculation  of  the  shear-work  integral  term,  Eq. 


(31). 


Figure  16.  Pressure  distributions  on  the  waisted 
body  of  revolution  for  = 0.597, 

Re.,  £ = 9.98  x 106  and  = 1.398, 

Re“'e  = 10.08  x 106  (Ref.  21). 

The  most  upstream  boundary-layer  profile  measurements  reported  in  Ref.  21  were  at 
x = 24  in.  Initial  conditions  near  the  nose  were  not  available.  This  is  usually  the  case  in 
practice;  therefore,  it  is  of  interest  to  investigate  the  influence  of  initial  conditions  on  the 
solutions.  Skin  friction  measurements  were  reported  for  x < 24  in.  in  Ref.  21  and  these 
measurements  provide  an  indication  of  the  results  that  might  be  expected.  Solutions  for 
four  sets  of  initial  conditions  are  compared  to  these  skin  friction  data  for  = 0.597 
and  Re^g  = 9.98  x 106  (Ref.  21)  in  Fig.  21.  Three  solutions  are  for  initial  conditions  at 
x = 3 in.  and  one  for  initial  conditions  at  x = 24  in.,  which  is  the  same  solution  as  given 
in  Fig.  17.  As  illustrated  in  Fig.  21,  the  three  solutions  beginning  at  x = 3 in.  quickly 
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x,  in. 

Figure  17.  Boundary-layer  properties  on  the  waisted 
body  of  revolution  for  = 0.597  and 
Re«,2  = 9 98  x 106. 

converge  to  the  same  axial  distribution  of  cf.  The  solution  beginning  at  x = 24  in.  also 
has  the  same  distribution,  and  converges  to  the  same  solution  as  the  three  solutions 
started  at  x = 3 in. 

The  data  of  Winter,  Rotta,  and  Smith  (Ref.  21)  were  used  for  comparisons  with 
other  calculation  techniques.  Blake  (Ref.  22)  compared  these  data  with  the  finite 
difference  calculation  methods  of  Cebeci  and  Smith  (Ref.  5),  Harris  (Ref.  23),  and  31ake 
(Ref.  22).  Three  figures  were  taken  from  Ref.  22  and  included  here  along  with  the 
present  calculations  as  Figs.  22,  23,  and  24.  It  is  important  to  note  that  these  figures 
were  copied  from  Ref.  22  because  the  plotted  experimental  data  points  differ  somewhat 
from  supposedly  the  same  data  (M^  = 1.398  and  Re^g  = 10.08  x 106) already  used  in 
Figs.  1 8 and  20.  The  discrepancies  are  not  large,  however,  and  direct  comparisons  of  the 
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calculation  methods  can  be  made.  Skin  friction  comparisons  in  Fig.  22  indicate  that  the 
present  solution  is  close  to  Blake's  solution  (Ref.  22)  up  to  about  x/2  = 0.7  and  then  is 
below  the  other  solutions  for  x/2  > 0.7  and  below  the  experimental  data  for  x/£  > 0.8. 
As  pointed  out  previously,  however,  it  was  concluded  by  Winter,  Rotta,  and  Smith  (Ref. 
21)  that  the  experimental  skin  friction  appears  too  large  over  the  aft  portion  of  the 
body.  Momentum  thickness  comparisons  in  Fig.  23  indicate  that  the  present  solution  is 
essentially  the  same  as  Cebeci  and  Smith  (Ref.  5)  up  to  the  maximum  value  of  0/Z,  and 
then  is  between  the  Cebeci  and  Smith  solution  and  Blake's  and  Harris'  solutions  (the  two 
latter  solutions  arc  essentially  the  same)  over  the  remaining  aft  end  of  the  body.  From 
the  velocity  distribution  comparisons  in  Fig.  24  the  present  solution  is  considered  to 
provide  a better  description  of  the  experimental  velocity  profile  than  Blake's  solution 
(Ref.  22).  Velocity  distributions  from  Cebeci  and  Smith  (Ref.  5)  and  Harris  (Ref.  23) 
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Figure  18.  Boundary-layer  properties  on  the 
waisted  body  of  revolution  for 
= 1.398  and  Re  { = 10.08 
x 10s. 
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calculation  methods  were  not  given  in  Ref.  22.  The  discrepancy  between  experimental 
data  and  the  present  solution  near  the  outer  edge  of  the  boundary  layer  was  discussed 
above  and  attributed  to  neglecting  transverse  curvature.  However,  the  slight  inflection  in 
the  experimental, velocity  profile  in  Fig.  24  is  well  represented  by  the  present  integral 
solution,  whereas  the  Finite  difference  solution  of  Blake  does  not  properly  describe  this 
feature  of  the  boundary  layer. 

The  experimental  data  of  Winter,  Rotta,  and  Smith  (Ref.  21)  were  also  considered 
by  Lewis,  Kubota,  and  Webb  (Ref.  9)  in  their  work  concerning  the  application  of 
transformation  theory  to  turbulent  boundary  layers  and  the  extension  of  Coles'  (Ref.  8) 
analysis  to  include  pressure  gradient.  The  agreement  obtained  in  Ref.  9 is  considerably 
worse  than  the  agreement  obtained  here  with  the  same  experimental  data.  Two  possible 
sources  of  discrepancy  are  given  in  Ref.  9 as  having  to  do  with  the  low-speed  formulation 
part  of  their  transformation  theory  and/or  the  Mangier  transformation  which  they  used. 
Such  problems  are  avoided  by  the  present  approach. 


Figure  19.  Velocity  profiles  on  the  waisted  body  of  revolution  for 
Wl  = 0.597  and  Re  £ = 9.98  x 106. 

oo  #0*4, 
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Figure  20.  Velocity  profiles  on  the  waisted  body  of  revolution  for 
M = 1.398  and  Re  £ = 10.08  x 106. 
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Figure  21.  Influence  of  initial  conditions  on  the 
waisted  body  of  revolution  for 

= 0.597  and  Re^  g = 9.98  x 106. 
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o Experimental  Data,  Ref.  21( 


Figure  22.  Calculated  and  measured  skin  friction 
distributions  on  the  waisted  body  of 
revolution  for  = 1 .398  and 
Re«,2  = 10  08  x°io«. 


Figure  23.  Calculated  and  measured  momentum 
thichness  distributions  on  the  waisted 
body  of  revolution  for  = 1.398 
arid  Re^g  = 10.08  x 106. 
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u/ue 

Figure  24.  Calculated  and  measured 
velocity  distributions  on 
the  waisted  body  of 
revolution  for  M„  = 1.398 
and  Re^g  = 10.08  x 106 
at  x = 42  in. 

4.0  CONCLUSIONS 

The  objective  of  the  present  work  was  the  development  of  a reasonably  accurate  and 
rapid  means  of  calculating  compressible  turbulent  adiabatic  wall  boundary  layers  on 
planar  and  axisymmetric  bodies.  Reasonably  accurate  solutions  were  obtained  for  the 
type  of  boundary-layer  problems  frequently  encountered  in  PWT  as  illustrated  by 
comparisons  with  transonic  turbulent  boundary-layer  data.  For  the  few  direct 
comparisons  afforded  by  the  results  in  Blake's  thesis  (Ref.  22),  the  present  integral 
method  provided  reasonable  agreement  with  the  finite  difference  methods  of  Cebeci  and 
Smith  (Ref.  5),  Blake  (Ref.  22),  and  Harris  (Ref.  23).  The  present  method  is  simple  to 
use  and  about  42  x-steps  can  be  solved  per  second  on  an  IBM  370/165.  The  method  no 
doubt  could  be  improved,  but  the  objective  of  the  present  work  is  considered  satisfied. 

The  present  approach  differs  from  most  integral  approaches  in  three  primary 
respects.  First,  an  entirely  different  velocity  profile  (Ref.  6)  was  used.  Second,  the 
shear-work  integral  term  was  not  correlated  from  experimental  data,  rather  it  was 
obtained  by  direct  integration.  And  third,  a new  relation  between  velocity  and 
temperature  was  used  for  shape  factor  calculations  (Ref.  14),  which  accounted  for 
total-temperature  overshoot  near  the  outer  edge  of  turbulent  boundary  layers  for  the 
realistic  conditions  of  nonunity  Prandtl  number  and  nonunity  recovery  factor. 
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APPENDIX  A 

LISTING  OF  THE  BOUNDARY-LAYER  COMPUTER  CODE 


IV  6 LEVEL  21  MAIN 

REAL  M.MINF,MM(300) .MOM (21 »KE(2) ,M2 
DIMENSION  X (2) . XX (30  0 1 .RR ( 300 > 

2 FORMAT  ( 1PUE12. 41 

3 FORMAT (8F 1 0 .5) 

4 FORMAT ( 13 ) 

6 FORMAT (///.3X, IMACH  NO.  ■ * . F6.3 . 5X . ■ UNI T REYNOLDS  NO.  *',F7.4,'  MI 
1LLI0N*.//) 

7 FORM AT(1H 1, 5X, IX* .1 IX. 'R*.11X.'M'.9X. 'THETA* .9X.*H*,11X.*CF*.9X. 

1 'DSTAR'.TX, ‘BETA* «9X. *REO* . 9X . * REX '.BX.'DISS'/I 

8 FORMAT!'  INPUT  VALUES  ARE • //6X . * XX • , 7X • * MM' ,8X . » RR « / 1 

9 FORMAT ( 3F 10 ,51 

12  READI5.31MINF, REN, THETA. H 
READ (5.4)  NINP 

READ (5.91 (XXII) ,MM(I) .RR(II .1=1, NINP) 

WRITE (6.8) 

WRITE (6.9) (XX (I) .MM(I) ,RR(I) .1=1, NINP) 

WRITE (6,6) M INF , REN 
X(])*XX<1) 

1*0 

IPR1NT*0 
500  KKK*0 

WRITE (6,7) 

1002  CONTINUE 

IF ( I ,LT .2) 1*1*1 
IF ( 1-1 >22,22,18 
18  DX*(XX(NINP)-XX(1) 1/500. 

20  X(I)*X(I-1) *DX 

IF ( (X ( I ) *DX-XX (NINP) ) .GE.0.0) IPRINT*1 
IF (X ( I )-XX (NINP) ) 22,22, 120 
22  CONTINUE 

CALL  LININT(XII) .XX.MM.M.DXM) 

CALL  LININTIX(I) ,XX,RR.RAD,DUM> 

M2*M*M 

T*(1.*0.2*MINF**2!/(1.»0,2*M**2) 

U*M/MINF*SQRT(T) 

DU*DXM/(M*(l.*0.2»M**2> ) 

FC*50RT(l.*M*M/5.) 

Q»U**2*T**2,5*RAD 

RE*REN*1.0E*06*U*T*M.T4 

HLI*1. 

IF(I-1)30. 30.36 
30  MOM(l)*Q*THETA 

HIb(H-0,290*M2)/(1.*o.113*M2) 

KE(1)*MOM(1)*U*(2.-0.54*TANH(1.1*(HI-1.) > *0.028*M2> / ( 1 ,*0.014*M2> 
REO*RE*MOM { 1 ) /Q/FC 

CF*0.3*EXP(-1.33*HI)/ALOG10(REO)** (1,74*0. 31 *HI) 

UPE*SQRT (2,/CF) 

CALL  TURB(M2,UPE,REQ.HI,THETA.IPRINT.O) 

36  DO  60  K*1 ■ 6 

IF (I“l 142,42,40 

40  MOM(I)*MOM(I-1)*0.5*(DM*OMM)«DX  . 

KE(H=KE(I-l)*0.5*(DK*DKK)*DX 
Z*(2.-KE(I)/MOM(I)/U)*(1.*0.014*M2)/0.54 
HI*l.*0.5*ALO6( (l.*Z)/(l.-Z) >/l.i 
H*HI*(1.*0.11 3*M2) *0.290*M2 

42  CF«0. 3*EXP (-1, 33*MI) /AL0G10 (RE«*MOM(I )/0/FC)** (1,74*0. 31*HI)/FC 
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FORTRAN 

0055 
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005S 
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0060 
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0065 
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0069 

0070 

0071 
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0078 
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0080 
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0084 


FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 
0009 


IV  e LEVEL  21  MAIN 

46  Dm«0.5*9*CF-M0M<I)*H*DU 

DK=Q*U*CF*D-2.*MOM<I)*U*DU*<0.064/<HI-0.8)*0»251)*M2 
IF  tl-l) 70»70»49 

49  IF (ABS<1. -HLI/HI 1-0.0005) 70 • 70 .50 

50  HLIbHI 
60  CONTINUE 

DX=0.5*DX 
GO  TO '20 
70  DMMbDM 
DKK"DK 

THETA«M0M(I)/0 
DSTAR  «H*THETA 
BETA»-2.*D$TAR/CF*DU 
REQCON»RE*THETA 
REX«RE*X<I> 

hRITE (6*2) X ( I ) t R AD *M«  THETA *HtCF*DSTARf BETA tREQCOM* REX tD 

upebsorti2./icf*fc)> 

REQbREQCOM/FC 

TMETAI«THETA/(l.-0.92*M2/(7.09*M214TANHn.49*CHI-0.9) I I 
CALL  TURB(M2,UPE.REQ.HI,THETAI,IPRINT.D> 

KKKbKKK»1 

I^(I.EO.l)  GO  TO  1002 
X(1)bX<2) 

HOHM)bMOM(2> 

KE(1)bKE(2) 

IFtKKK.EQ. 57)00  TO  500 
GO  TO  18 
120  CONTINUE 
STOP 
END 


IV  G LEVEL  21  LIN1NT 

SUBROUTINE  LlNINT (X* XX, YY ,Y, DY» 

DIMENSION  XX (300) • YY (300 1 

Nsl 

1 N«N* 1 

IF(X.GT.XXIN))  GO  TO  1 

DYs(YY(N)-YY(N-l) )/(XX(N)-XX(N-l>> 

Y«YY (N-l ) * (X-XX (N- 1 ) )*DY 

RETURN 

END 
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0040 
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0044 


IV  G LEVEL  21  TURB 

SUBROUTINE  TURBIXM2.UPE. RET, H, THETA, IPRINT.O) 

DIMENSION  F (101) 

Pl"3. 141592 
RETUP=RET/(UPE**2> 

HEUP9=0.09*RET/UPE 

PIUPE»PI/(.18«UPE) 

UUE2=1.723*EXP(-0.6*H)*(1.*50./HET) 

UUE5«0.87*0.08*EXP(-2.6*IH-1.95)**2> 

G2« (UUE2* ( 1 ./ ( ,09*UPE> ) *ATAN (2,*REUP9) ) / ( 1 ,-PIUPEJ 
G5- (UUE5- ( 1 ./ t ,09*UPE> > *ATAN  <5.*REUP9) ) t < 1.-P1UPE) 
BaAL0G((,5*AL0G(< 1 «*G2**2) / < 1 ,-G2**2) > ) / < ( ,5*AL0G ( ( 1 •♦85**2) X II •" 
165**2) )) >) /ALOGIO. 4) 

A»(.5*AL0G(<l.*G2**2)/ll.-G2**2) ) 1X2. **B 
N»101 

DXelO./lOO. 

SUM.O. 

X«UPE*100./RET-DX 

DO  9000  I "1 *N 

X«X*DX- 

V»X*THETA 

VPLUS«X*RET/UPE 

ARG«A*X**B 

IF(ARG.GT.8T.3)F(I>»0. 

IF(ARG.GT.87.3)G0  TO  9000 
NCOS® (COSH ( ARG) ) **2 
HTANaSQRT (TANH(ARG) > 

UPLUS»ATAN(REUP9»X)/0.09*  (UP.E-PI/0. 1 8)  *HTAN 
UUE«UPLUS/UPE 

F ( I ) "RETUP* Il./( 1.* |REUP9*X) **2) ) ♦ < 1 ,-PIUPE) *ARG*B/2./X 
1 * < 1 • / ( HT AN*HC0S) ) 

XL-0.41*X 

FFC«0.016B*H*F(I) 

FF«XL*XL*F(I>*F (1) 

IF (FF.GT.FFC) SUM=SUM*1  . 

IFISUM.GT.O.)FF»FFC 

F ( I ) «FF*F ( I ) / ( 1 . *0 . 1 76*XM2« 1 1 .-UUE*UUE ) ) 
IF(IPRINT.EQ.l.AND.I.EQ.l)  WRITE (6 , 10) 

10  FORMAT (,l,*5X**Y*i  8X. « V/THETAI * , 7X, • Y* • « lOXi *U* ' • 9X,«U/UE*/) 

IF < IPRINT «EQ, 1 ) WRITE(6<11)  V.X. YPLUS.UPLUS.UUE 

11  FORMAT ( 1P5E12.4) 

9000  CONTINUE 

CALL  SIMP IF,F2,F4,N) 

CP»DX/3.*(F(1) *F2*F4*F  IN) ) 

D=16.79/UPE*UPE*UPE*CP 

RETURN 

END 
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FORTRAN  IV  G LEVEL  21  SIMP 


0001 

SUBROUTINE  SIMP(F,F2.F4.N) 

0002 

DIMENSION  F (101) 

0003 

FODDoO.O 

0004 

FEVENbO.O 

0005 

N1«N-1 

0006 

N2-N-2 

0007 

DO  30  I»3tN2«2 

ooos 

30 

FODDoFODO^F (I) 

0009 

F2-2.«FODD 

0010 

DO  40  I*2«N1 12 

0011 

40 

FEVEN«FEVEN*F ( I ) 

0012 

F4»4.*FEVEN 

0013 

RETURN 

0014 

END 

The  following  is  an  example  set  of  input  data  for  the  Waisted  Body  of  Revolution  (Ref. 
21)  for  = 0.597  and  Re„,g  = 9.98  x 10<\ 


.597 

022 

.1642 

.0363 

24.0 

.642 

4.358 

27.0 

.598 

3.793 

2B.5 

.585 

3.515 

30.0 

.578 

3.257 

31.5 

.573 

3.0199 

33.0 

.569 

2. 7981 

34.5 

.567 

2.544 

36.0 

.563  • 

2.404 

37.5 

.559 

2.222 

40.0 

.545 

1.969 

42.0 

.537 

1.87 

44.0 

.542 

1.933 

45.5 

.552 

2.067 

47.0 

.566 

2.231 

48.5 

.663 

2.407 

50.0 

.567 

2.597 

51.5 

.569 

2.80? 

53.0 

.573 

3.023 

54.5 

.578 

3.262 

56.0 

.585 

1.52 

57.5 

.597 

3.79? 

59.0 

.623 

4.097 
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NOMENCLATURE 

a Parameter  used  and  defined  in  Table  1 

b Parameter  used  and  defined  in  Table  1 

c Airfoil  chord 

Cf  Local  skin  friction  coefficient,  /qe 

cp  Pressure  coefficient,  (p-pj/q^ 

cp*  Critical  pressure  coefficient 

D Shear-work  integral,  defined  by  Eq.  (26) 

D Defined  by  Eq.  (19) 

Fc  (1  + 0.2  Mg)1 

g Function  used  and  defined  in  Table  1 

H Specific  total  enthalpy 

H 8*16 

Hg*  5*/0 

h8**  5 **ie 

Hg*  e*/e 

k Index  which  is  zero  for  planar  flow  and  one  for  axisymmetric  flow 

S.  Mixing  length,  £ = 0.41  y;  also  total  body  length 

M Mach  number 

m Index  which  is  zero  for  the  momentum  integral  equation  and  one  for  the 

mean-flow  kinetic  energy  integral  equation 

Pr  Molecular  Prandtl  number 

Prt  Turbulent  Prandtl  number 
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p Static  pressure 

q Dynamic  pressure.  pu2/2 

Reg  pL.  i\.  0/pL. 

Re„,L  p„  u_  c/p„ 

Re-.e  p.  u„ 

r rw.  ± y cos  \js;  also  used  as  recovery  factor  (r  = 0.88) 

rw  Local  body  radius 

T Static  temperature 

u Mean  velocity  in  the  x direction 

ll+  u/uT 

uT  (rw/Pw)1/2 

v Mean  velocity  in  the  y direction 

x Coordinate  along  the  body  surface 

y Coordinate  normal  to  the  body  surface 

- k. 

y+  Pw  uT  y/pw . (Re0/ue+)(y !6) 

a Angle  of  attack,  deg 

(3  5*/tw  dp/dx 

5*  Boundary-layer  displacement  thickness  defined  by  Eq.  (23) 

8*  Defined  by  Eq.  (36) 

8*  Defined  by  Eq.  (16) 

5**  Defined  by  Eq.  (25) 

5**  Defined  by  Eq.  (18) 

e Eddy  viscosity,  defined  by  Eq.  (32)  i 

em  Middle  boundary-layer  region  eddy  viscosity,  defined  by  Eq.  (33) 
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eo  Outer  boundary-layer  region  eddy  viseosity,  defined  by  Eq.  (35) 

8 Boundary-layer  momentum  thickness,  defined  by  Eq.  (22) 

6 Low-speed  or  kinematic  boundary-layer  momentum  thickness,  B = /“u/ue 

(l-u/ue)dy  0 

6 Defined  by  Eq.  (15) 

8*  Defined  by  Eq.  (24) 

8*  Defined  by  Eq.  (17) 

fi  Molecular  viscosity 

p Density 

-pu'v1  Reynolds  stress 

t Total  shear  stress,  (p  + pe)  3u/3y 

rt  Turbulent  shear  stress,  pe  3u/dy 

i//  Local  body  surface  angle 

SUBSCRIPTS 

e Boundary-layer  edge  value 

exp  Experimental  value 

w Wall  value 

Free-stream  value 

SUPERSCRIPT 

— Denotes  low-speed  or  incompressible  value 
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